library(tidyverse) # tidy universe
library(patchwork) # combine plots
library(lmerTest) # mixed model
library(car) # ANOVA
library(emmeans) # post hoc
library(performance) # model performanceConfounder boar (Igf & animal welfare)
R Libraries
set.seed(1989)Data
Read data
Data processing done in file “igf-biomarker-summary-stats.qmd”.
load("./data/data-processed.RData")table(dat.w$husbandry, dat.w$boar)
Astein Bestein Beton Capello Faro Flausist Sacher Slavon Slovan
conventional 1 0 6 6 0 2 2 1 2
ecological 4 2 9 1 1 1 0 0 0
Tonkar
conventional 0
ecological 2Statistical model: general design
With a simple model structure we solely test for a confounder effect of the boars. We include a random animal effect and fit a linear mixed model with the lmerTest package in R.
contr = lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 10000000),
calc.derivs = FALSE)Birth weight
Model
hist(dat.w$mean.bodyweight,
breaks = 30)
mod.bodyweight <- lmerTest::lmer(mean.bodyweight ~
boar +
# random intercept for sows
(1 | sow),
data = dat.w,
REML = TRUE,
control = contr)summary(mod.bodyweight)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: mean.bodyweight ~ boar + (1 | sow)
Data: dat.w
Control: contr
REML criterion at convergence: 10.2
Scaled residuals:
Min 1Q Median 3Q Max
-2.23173 -0.57132 0.08495 0.58636 1.93133
Random effects:
Groups Name Variance Std.Dev.
sow (Intercept) 0.003039 0.05513
Residual 0.055767 0.23615
Number of obs: 40, groups: sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.28340 0.10830 29.99958 11.850 7.63e-13 ***
boarBestein 0.19258 0.20261 29.99985 0.951 0.349
boarBeton 0.12817 0.12416 27.88475 1.032 0.311
boarCapello -0.05304 0.14080 27.95471 -0.377 0.709
boarFaro -0.16233 0.26503 29.87555 -0.612 0.545
boarFlausist 0.04580 0.17568 28.20866 0.261 0.796
boarSacher 0.28907 0.20241 29.86622 1.428 0.164
boarSlavon 0.08854 0.26514 29.96194 0.334 0.741
boarSlovan -0.11074 0.20261 29.99998 -0.547 0.589
boarTonkar 0.04312 0.19963 23.10452 0.216 0.831
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) brBstn borBtn brCpll boarFr brFlss brSchr boarSlavn
boarBestein -0.534
boarBeton -0.860 0.467
boarCapello -0.752 0.402 0.654
boarFaro -0.408 0.218 0.363 0.307
boarFlausst -0.606 0.324 0.527 0.462 0.248
boarSacher -0.534 0.286 0.467 0.417 0.218 0.324
boarSlavon -0.408 0.218 0.357 0.307 0.167 0.248 0.218
boarSlovan -0.534 0.285 0.464 0.402 0.218 0.338 0.285 0.245
boarTonkar -0.515 0.276 0.451 0.388 0.211 0.312 0.276 0.211
boarSlovn
boarBestein
boarBeton
boarCapello
boarFaro
boarFlausst
boarSacher
boarSlavon
boarSlovan
boarTonkar 0.275 round(drop1(mod.bodyweight, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
mean.bodyweight ~ boar + (1 | sow)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
boar 0.409 0.045 9 28.805 0.815 0.606car::Anova(mod.bodyweight,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: mean.bodyweight
Chisq Df Pr(>Chisq)
boar 7.3391 9 0.6019car::Anova(mod.bodyweight,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: mean.bodyweight
F Df Df.res Pr(>F)
boar 0.7149 9 26.965 0.6907Model diagnostics
performance::check_model(mod.bodyweight)
Plot
plot.bodyweight <- dat.w |>
# make plot
mutate(jit = jitter(as.numeric(boar), 0.3)) |>
ggplot(aes(y = mean.bodyweight)) +
geom_boxplot(aes(x = boar,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = boar,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 2),
breaks = seq(0, 2, 0.5),
minor_breaks = seq(0, 2, by = 0.1) ) +
labs(x = "",
y = "Mean birth weight [kg]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.bodyweight
Total piglets
Model
hist(dat.w$total.piglets,
breaks = 30)
mod.piglets <- lmerTest::lmer(total.piglets ~
boar +
# random intercept for sows
(1 | sow),
data = dat.w,
REML = TRUE,
control = contr)summary(mod.piglets)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: total.piglets ~ boar + (1 | sow)
Data: dat.w
Control: contr
REML criterion at convergence: 173
Scaled residuals:
Min 1Q Median 3Q Max
-2.2833 -0.4447 0.0667 0.5248 1.3179
Random effects:
Groups Name Variance Std.Dev.
sow (Intercept) 5.536 2.353
Residual 9.320 3.053
Number of obs: 40, groups: sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 17.3468 1.6209 29.6138 10.702 1.08e-11 ***
boarBestein -3.4548 3.0186 28.3752 -1.145 0.262
boarBeton -0.2354 1.7414 16.6375 -0.135 0.894
boarCapello -1.2092 1.9754 16.5531 -0.612 0.549
boarFaro 5.3426 3.7983 18.4212 1.407 0.176
boarFlausist -1.8109 2.4825 18.2299 -0.729 0.475
boarSacher -1.5936 2.8986 18.2757 -0.550 0.589
boarSlavon 3.4623 3.8503 21.4644 0.899 0.379
boarSlovan 3.2054 2.9944 25.6667 1.070 0.294
boarTonkar 0.5989 2.6953 12.5549 0.222 0.828
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) brBstn borBtn brCpll boarFr brFlss brSchr boarSlavn
boarBestein -0.516
boarBeton -0.809 0.477
boarCapello -0.676 0.361 0.616
boarFaro -0.396 0.219 0.407 0.285
boarFlausst -0.562 0.300 0.511 0.431 0.237
boarSacher -0.513 0.276 0.484 0.458 0.220 0.307
boarSlavon -0.400 0.216 0.376 0.282 0.172 0.256 0.216
boarSlovan -0.514 0.273 0.462 0.363 0.215 0.401 0.274 0.367
boarTonkar -0.417 0.226 0.398 0.294 0.180 0.244 0.226 0.176
boarSlovn
boarBestein
boarBeton
boarCapello
boarFaro
boarFlausst
boarSacher
boarSlavon
boarSlovan
boarTonkar 0.222 round(drop1(mod.piglets, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
total.piglets ~ boar + (1 | sow)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
boar 78.053 8.673 9 18.337 0.93 0.522car::Anova(mod.piglets,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: total.piglets
Chisq Df Pr(>Chisq)
boar 8.3744 9 0.4969car::Anova(mod.piglets,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: total.piglets
F Df Df.res Pr(>F)
boar 0.84 9 23.759 0.5877Model diagnostics
performance::check_model(mod.piglets)
Plot
plot.piglets <- dat.w |>
# make plot
mutate(jit = jitter(as.numeric(boar), 0.3)) |>
ggplot(aes(y = total.piglets)) +
geom_boxplot(aes(x = boar,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = boar,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 22),
breaks = seq(0, 22, 5),
minor_breaks = seq(0, 22, by = 1) ) +
labs(x = "",
y = "Number of piglets",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.piglets
Proportion males
Model
hist(dat.w$prop.males,
breaks = 30)
mod.males <- lmerTest::lmer(prop.males ~
boar +
# random intercept for sows
(1 | sow),
data = dat.w,
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.males)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: prop.males ~ boar + (1 | sow)
Data: dat.w
Control: contr
REML criterion at convergence: 246.1
Scaled residuals:
Min 1Q Median 3Q Max
-2.3908 -0.4465 0.0000 0.5513 2.2803
Random effects:
Groups Name Variance Std.Dev.
sow (Intercept) 0.0 0.00
Residual 152.7 12.36
Number of obs: 40, groups: sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 55.3165 5.5268 30.0000 10.009 4.48e-11 ***
boarBestein 0.4527 10.3397 30.0000 0.044 0.965
boarBeton -10.7701 6.3818 30.0000 -1.688 0.102
boarCapello -8.9729 7.2363 30.0000 -1.240 0.225
boarFaro -17.2213 13.5379 30.0000 -1.272 0.213
boarFlausist -1.2389 9.0253 30.0000 -0.137 0.892
boarSacher -10.0959 10.3397 30.0000 -0.976 0.337
boarSlavon 2.5782 13.5379 30.0000 0.190 0.850
boarSlovan 0.9993 10.3397 30.0000 0.097 0.924
boarTonkar 1.5548 10.3397 30.0000 0.150 0.881
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) brBstn borBtn brCpll boarFr brFlss brSchr boarSlavn
boarBestein -0.535
boarBeton -0.866 0.463
boarCapello -0.764 0.408 0.661
boarFaro -0.408 0.218 0.354 0.312
boarFlausst -0.612 0.327 0.530 0.468 0.250
boarSacher -0.535 0.286 0.463 0.408 0.218 0.327
boarSlavon -0.408 0.218 0.354 0.312 0.167 0.250 0.218
boarSlovan -0.535 0.286 0.463 0.408 0.218 0.327 0.286 0.218
boarTonkar -0.535 0.286 0.463 0.408 0.218 0.327 0.286 0.218
boarSlovn
boarBestein
boarBeton
boarCapello
boarFaro
boarFlausst
boarSacher
boarSlavon
boarSlovan
boarTonkar 0.286
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.males, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
prop.males ~ boar + (1 | sow)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
boar 1172.5 130.28 9 30 0.853 0.575car::Anova(mod.males,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: prop.males
Chisq Df Pr(>Chisq)
boar 7.677 9 0.567car::Anova(mod.males,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: prop.males
F Df Df.res Pr(>F)
boar 0.733 9 26.984 0.6756Model diagnostics
performance::check_model(mod.males)
Plot
plot.males <- dat.w |>
# make plot
mutate(jit = jitter(as.numeric(boar), 0.3)) |>
ggplot(aes(y = prop.males)) +
geom_boxplot(aes(x = boar,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = boar,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 80),
breaks = seq(0, 80, 20),
minor_breaks = seq(0, 80, by = 5) ) +
labs(x = "",
y = "Proportion males [%]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.males
Proportion stillborn
Model
hist(dat.w$prop.stillborn,
breaks = 30)
hist(log(dat.w$prop.stillborn+1),
breaks = 30)
mod.stillborn <- lmerTest::lmer(log(dat.w$prop.stillborn+1) ~
boar +
# random intercept for sows
(1 | sow),
data = dat.w,
REML = TRUE,
control = contr)summary(mod.stillborn)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(dat.w$prop.stillborn + 1) ~ boar + (1 | sow)
Data: dat.w
Control: contr
REML criterion at convergence: 106.9
Scaled residuals:
Min 1Q Median 3Q Max
-1.32881 -0.42267 0.00086 0.47963 1.82883
Random effects:
Groups Name Variance Std.Dev.
sow (Intercept) 0.426 0.6527
Residual 1.138 1.0668
Number of obs: 40, groups: sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.25893 0.54075 29.86540 2.328 0.0269 *
boarBestein -1.11985 1.01007 29.69135 -1.109 0.2765
boarBeton 0.31743 0.59564 24.42302 0.533 0.5989
boarCapello -0.06806 0.67573 24.42948 -0.101 0.9206
boarFaro 1.72993 1.29344 26.18582 1.337 0.1926
boarFlausist 1.25723 0.84656 25.36699 1.485 0.1498
boarSacher -0.20312 0.98734 26.09434 -0.206 0.8386
boarSlavon -0.79010 1.30408 27.67283 -0.606 0.5495
boarSlovan 1.14679 1.00678 29.18801 1.139 0.2639
boarTonkar -1.30818 0.93170 20.48662 -1.404 0.1753
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) brBstn borBtn brCpll boarFr brFlss brSchr boarSlavn
boarBestein -0.525
boarBeton -0.828 0.477
boarCapello -0.702 0.374 0.626
boarFaro -0.402 0.219 0.395 0.291
boarFlausst -0.577 0.308 0.516 0.440 0.240
boarSacher -0.523 0.281 0.480 0.447 0.220 0.312
boarSlavon -0.405 0.217 0.371 0.290 0.170 0.250 0.217
boarSlovan -0.524 0.279 0.463 0.375 0.217 0.382 0.279 0.335
boarTonkar -0.445 0.239 0.412 0.319 0.188 0.262 0.239 0.185
boarSlovn
boarBestein
boarBeton
boarCapello
boarFaro
boarFlausst
boarSacher
boarSlavon
boarSlovan
boarTonkar 0.237 round(drop1(mod.stillborn, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(dat.w$prop.stillborn + 1) ~ boar + (1 | sow)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
boar 14.102 1.567 9 25.813 1.377 0.249car::Anova(mod.stillborn,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(dat.w$prop.stillborn + 1)
Chisq Df Pr(>Chisq)
boar 12.391 9 0.1922car::Anova(mod.stillborn,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(dat.w$prop.stillborn + 1)
F Df Df.res Pr(>F)
boar 1.2425 9 25.027 0.3149Model diagnostics
performance::check_model(mod.stillborn)
Plot
plot.stillborn <- dat.w |>
# make plot
mutate(jit = jitter(as.numeric(boar), 0.3)) |>
ggplot(aes(y = prop.stillborn)) +
geom_boxplot(aes(x = boar,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = boar,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 40),
breaks = seq(0, 40, 10),
minor_breaks = seq(0, 40, by = 2) ) +
labs(x = "",
y = "Proportion stillborn [%]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.stillborn
Proportion spreizer
Model
hist(dat.w$prop.spreizer,
breaks = 30)
hist(log(dat.w$prop.spreizer+1),
breaks = 30)
mod.spreizer <- lmerTest::lmer(log(dat.w$prop.spreizer+1) ~
boar +
# random intercept for sows
(1 | sow),
data = dat.w,
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.spreizer)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(dat.w$prop.spreizer + 1) ~ boar + (1 | sow)
Data: dat.w
Control: contr
REML criterion at convergence: 79.6
Scaled residuals:
Min 1Q Median 3Q Max
-1.3649 -0.1482 -0.1482 0.0000 2.4683
Random effects:
Groups Name Variance Std.Dev.
sow (Intercept) 0.0000 0.0000
Residual 0.5938 0.7706
Number of obs: 40, groups: sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 9.692e-16 3.446e-01 3.000e+01 0.000 1.000000
boarBestein -1.128e-15 6.447e-01 3.000e+01 0.000 1.000000
boarBeton 1.142e-01 3.979e-01 3.000e+01 0.287 0.776102
boarCapello 1.052e+00 4.512e-01 3.000e+01 2.331 0.026660 *
boarFaro 1.751e+00 8.441e-01 3.000e+01 2.075 0.046699 *
boarFlausist 6.603e-01 5.628e-01 3.000e+01 1.173 0.249878
boarSacher -9.460e-16 6.447e-01 3.000e+01 0.000 1.000000
boarSlavon -9.362e-16 8.441e-01 3.000e+01 0.000 1.000000
boarSlovan 2.609e+00 6.447e-01 3.000e+01 4.046 0.000336 ***
boarTonkar 9.173e-01 6.447e-01 3.000e+01 1.423 0.165100
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) brBstn borBtn brCpll boarFr brFlss brSchr boarSlavn
boarBestein -0.535
boarBeton -0.866 0.463
boarCapello -0.764 0.408 0.661
boarFaro -0.408 0.218 0.354 0.312
boarFlausst -0.612 0.327 0.530 0.468 0.250
boarSacher -0.535 0.286 0.463 0.408 0.218 0.327
boarSlavon -0.408 0.218 0.354 0.312 0.167 0.250 0.218
boarSlovan -0.535 0.286 0.463 0.408 0.218 0.327 0.286 0.218
boarTonkar -0.535 0.286 0.463 0.408 0.218 0.327 0.286 0.218
boarSlovn
boarBestein
boarBeton
boarCapello
boarFaro
boarFlausst
boarSacher
boarSlavon
boarSlovan
boarTonkar 0.286
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.spreizer, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(dat.w$prop.spreizer + 1) ~ boar + (1 | sow)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
boar 17.747 1.972 9 30 3.321 0.006 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.spreizer,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(dat.w$prop.spreizer + 1)
Chisq Df Pr(>Chisq)
boar 29.886 9 0.0004587 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.spreizer,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(dat.w$prop.spreizer + 1)
F Df Df.res Pr(>F)
boar 3.0219 9 26.984 0.01249 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model diagnostics
performance::check_model(mod.spreizer)
Emmeans & Effect sizes
Emmeans:
emm <- emmeans(mod.spreizer,
pairwise ~ boar,
data = dat.w,
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
boar response SE df lower.CL upper.CL
Astein 0.000 0.345 30 -0.505 1.021
Bestein 0.000 0.545 30 -0.671 2.043
Beton 0.121 0.223 30 -0.253 0.683
Capello 1.863 0.834 30 0.579 4.189
Faro 4.762 4.440 30 0.194 26.800
Flausist 0.935 0.861 30 -0.220 3.802
Sacher 0.000 0.545 30 -0.671 2.043
Slavon 0.000 0.771 30 -0.793 3.825
Slovan 12.580 7.400 30 3.463 40.324
Tonkar 1.503 1.360 30 -0.178 6.615
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log(mu + 1) scale
$contrasts
contrast estimate SE df t.ratio p.value
Astein - Bestein 0.000 0.645 30 0.000 1.0000
Astein - Beton -0.114 0.398 30 -0.287 1.0000
Astein - Capello -1.052 0.451 30 -2.331 0.3995
Astein - Faro -1.751 0.844 30 -2.075 0.5576
Astein - Flausist -0.660 0.563 30 -1.173 0.9711
Astein - Sacher 0.000 0.645 30 0.000 1.0000
Astein - Slavon 0.000 0.844 30 0.000 1.0000
Astein - Slovan -2.609 0.645 30 -4.046 0.0106
Astein - Tonkar -0.917 0.645 30 -1.423 0.9103
Bestein - Beton -0.114 0.580 30 -0.197 1.0000
Bestein - Capello -1.052 0.618 30 -1.702 0.7856
Bestein - Faro -1.751 0.944 30 -1.856 0.6966
Bestein - Flausist -0.660 0.703 30 -0.939 0.9936
Bestein - Sacher 0.000 0.771 30 0.000 1.0000
Bestein - Slavon 0.000 0.944 30 0.000 1.0000
Bestein - Slovan -2.609 0.771 30 -3.385 0.0531
Bestein - Tonkar -0.917 0.771 30 -1.190 0.9683
Beton - Capello -0.938 0.353 30 -2.658 0.2357
Beton - Faro -1.637 0.796 30 -2.057 0.5689
Beton - Flausist -0.546 0.487 30 -1.121 0.9785
Beton - Sacher 0.114 0.580 30 0.197 1.0000
Beton - Slavon 0.114 0.796 30 0.143 1.0000
Beton - Slovan -2.494 0.580 30 -4.300 0.0055
Beton - Tonkar -0.803 0.580 30 -1.385 0.9227
Capello - Faro -0.699 0.824 30 -0.849 0.9969
Capello - Flausist 0.391 0.532 30 0.736 0.9990
Capello - Sacher 1.052 0.618 30 1.702 0.7856
Capello - Slavon 1.052 0.824 30 1.277 0.9514
Capello - Slovan -1.557 0.618 30 -2.520 0.2985
Capello - Tonkar 0.134 0.618 30 0.218 1.0000
Faro - Flausist 1.091 0.890 30 1.226 0.9619
Faro - Sacher 1.751 0.944 30 1.856 0.6966
Faro - Slavon 1.751 1.090 30 1.607 0.8345
Faro - Slovan -0.857 0.944 30 -0.908 0.9950
Faro - Tonkar 0.834 0.944 30 0.884 0.9959
Flausist - Sacher 0.660 0.703 30 0.939 0.9936
Flausist - Slavon 0.660 0.890 30 0.742 0.9989
Flausist - Slovan -1.948 0.703 30 -2.770 0.1924
Flausist - Tonkar -0.257 0.703 30 -0.365 1.0000
Sacher - Slavon 0.000 0.944 30 0.000 1.0000
Sacher - Slovan -2.609 0.771 30 -3.385 0.0531
Sacher - Tonkar -0.917 0.771 30 -1.190 0.9683
Slavon - Slovan -2.609 0.944 30 -2.764 0.1944
Slavon - Tonkar -0.917 0.944 30 -0.972 0.9918
Slovan - Tonkar 1.691 0.771 30 2.195 0.4814
Note: contrasts are still on the log(mu + 1) scale. Consider using
regrid() if you want contrasts of back-transformed estimates.
Degrees-of-freedom method: satterthwaite
P value adjustment: tukey method for comparing a family of 10 estimates Effect sizes:
eff_size(emm,
sigma = sigma(mod.spreizer),
edf = df.residual(mod.spreizer))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(Astein - Bestein) 0.000 0.837 30 -1.7087 1.7087
(Astein - Beton) -0.148 0.517 30 -1.2036 0.9072
(Astein - Capello) -1.365 0.613 30 -2.6174 -0.1124
(Astein - Faro) -2.273 1.140 30 -4.5942 0.0490
(Astein - Flausist) -0.857 0.739 30 -2.3666 0.6528
(Astein - Sacher) 0.000 0.854 30 -1.7450 1.7450
(Astein - Slavon) 0.000 1.110 30 -2.2651 2.2651
(Astein - Slovan) -3.385 0.951 30 -5.3276 -1.4427
(Astein - Tonkar) -1.190 0.852 30 -2.9297 0.5489
(Bestein - Beton) -0.148 0.753 30 -1.6861 1.3897
(Bestein - Capello) -1.365 0.822 30 -3.0442 0.3144
(Bestein - Faro) -2.273 1.260 30 -4.8496 0.3044
(Bestein - Flausist) -0.857 0.920 30 -2.7359 1.0220
(Bestein - Sacher) 0.000 1.010 30 -2.0728 2.0728
(Bestein - Slavon) 0.000 1.240 30 -2.5262 2.5262
(Bestein - Slovan) -3.385 1.100 30 -5.6267 -1.1437
(Bestein - Tonkar) -1.190 1.010 30 -3.2584 0.8775
(Beton - Capello) -1.217 0.486 30 -2.2088 -0.2247
(Beton - Faro) -2.124 1.070 30 -4.3119 0.0631
(Beton - Flausist) -0.709 0.640 30 -2.0148 0.5973
(Beton - Sacher) 0.148 0.753 30 -1.3897 1.6861
(Beton - Slavon) 0.148 1.030 30 -1.9614 2.2578
(Beton - Slovan) -3.237 0.868 30 -5.0101 -1.4639
(Beton - Tonkar) -1.042 0.766 30 -2.6057 0.5212
(Capello - Faro) -0.908 1.080 30 -3.1050 1.2896
(Capello - Flausist) 0.508 0.693 30 -0.9081 1.9241
(Capello - Sacher) 1.365 0.822 30 -0.3144 3.0442
(Capello - Slavon) 1.365 1.080 30 -0.8499 3.5798
(Capello - Slovan) -2.020 0.846 30 -3.7480 -0.2925
(Capello - Tonkar) 0.174 0.802 30 -1.4637 1.8126
(Faro - Flausist) 1.416 1.170 30 -0.9740 3.8054
(Faro - Sacher) 2.273 1.260 30 -0.3044 4.8496
(Faro - Slavon) 2.273 1.450 30 -0.6814 5.2267
(Faro - Slovan) -1.113 1.230 30 -3.6322 1.4071
(Faro - Tonkar) 1.082 1.230 30 -1.4365 3.6008
(Flausist - Sacher) 0.857 0.920 30 -1.0220 2.7359
(Flausist - Slavon) 0.857 1.160 30 -1.5129 3.2267
(Flausist - Slovan) -2.528 0.973 30 -4.5162 -0.5403
(Flausist - Tonkar) -0.334 0.914 30 -2.2001 1.5330
(Sacher - Slavon) 0.000 1.240 30 -2.5262 2.5262
(Sacher - Slovan) -3.385 1.100 30 -5.6267 -1.1437
(Sacher - Tonkar) -1.190 1.010 30 -3.2584 0.8775
(Slavon - Slovan) -3.385 1.310 30 -6.0516 -0.7188
(Slavon - Tonkar) -1.190 1.240 30 -3.7127 1.3318
(Slovan - Tonkar) 2.195 1.040 30 0.0665 4.3230
sigma used for effect sizes: 0.7706
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Plot
plot.spreizer <- dat.w |>
# make plot
mutate(jit = jitter(as.numeric(boar), 0.3)) |>
ggplot(aes(y = prop.spreizer)) +
geom_boxplot(aes(x = boar,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = boar,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 20),
breaks = seq(0, 20, 5),
minor_breaks = seq(0, 20, by = 1) ) +
labs(x = "",
y = "Proportion spreizer [%]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.spreizer
Proportion low body weight
Model
hist(dat.w$prop.lbw,
breaks = 30)
hist(log(dat.w$prop.lbw+1),
breaks = 30)
mod.lbw <- lmerTest::lmer(log(dat.w$prop.lbw+1) ~
boar +
# random intercept for sows
(1 | sow),
data = dat.w,
REML = TRUE,
control = contr)summary(mod.lbw)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(dat.w$prop.lbw + 1) ~ boar + (1 | sow)
Data: dat.w
Control: contr
REML criterion at convergence: 112.5
Scaled residuals:
Min 1Q Median 3Q Max
-1.86969 -0.18804 0.02804 0.51646 1.57564
Random effects:
Groups Name Variance Std.Dev.
sow (Intercept) 0.03886 0.1971
Residual 1.73722 1.3180
Number of obs: 40, groups: sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 3.2141 0.5959 29.9846 5.394 7.67e-06 ***
boarBestein -1.0909 1.1147 29.9854 -0.979 0.336
boarBeton -0.7084 0.6860 28.0416 -1.033 0.311
boarCapello 0.2217 0.7779 28.1133 0.285 0.778
boarFaro 0.6673 1.4593 29.9971 0.457 0.651
boarFlausist -1.1866 0.9703 28.3199 -1.223 0.231
boarSacher -1.2459 1.1145 29.9964 -1.118 0.272
boarSlavon -0.3499 1.4594 29.9994 -0.240 0.812
boarSlovan 0.4711 1.1147 29.9849 0.423 0.676
boarTonkar -0.2093 1.1078 22.6429 -0.189 0.852
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) brBstn borBtn brCpll boarFr brFlss brSchr boarSlavn
boarBestein -0.534
boarBeton -0.863 0.465
boarCapello -0.759 0.406 0.658
boarFaro -0.408 0.218 0.358 0.310
boarFlausst -0.610 0.326 0.529 0.465 0.249
boarSacher -0.534 0.286 0.465 0.412 0.218 0.326
boarSlavon -0.408 0.218 0.355 0.310 0.167 0.249 0.218
boarSlovan -0.534 0.286 0.463 0.406 0.218 0.332 0.286 0.230
boarTonkar -0.526 0.281 0.458 0.399 0.215 0.321 0.281 0.215
boarSlovn
boarBestein
boarBeton
boarCapello
boarFaro
boarFlausst
boarSacher
boarSlavon
boarSlovan
boarTonkar 0.281 round(drop1(mod.lbw, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(dat.w$prop.lbw + 1) ~ boar + (1 | sow)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
boar 11.706 1.301 9 28.663 0.749 0.662car::Anova(mod.lbw,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(dat.w$prop.lbw + 1)
Chisq Df Pr(>Chisq)
boar 6.7386 9 0.6643car::Anova(mod.lbw,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(dat.w$prop.lbw + 1)
F Df Df.res Pr(>F)
boar 0.6549 9 27.009 0.7409Model diagnostics
performance::check_model(mod.lbw)
Plot
plot.lbw <- dat.w |>
# make plot
mutate(jit = jitter(as.numeric(boar), 0.3)) |>
ggplot(aes(y = prop.lbw)) +
geom_boxplot(aes(x = boar,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = boar,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 82),
breaks = seq(0, 82, 20),
minor_breaks = seq(0, 82, by = 5) ) +
labs(x = "",
y = "Proportion LBW [%]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.lbw
Combined plots: Figure
# Combine plots with a designated area for the legend
combined <- (plot.bodyweight +
plot.piglets +
guide_area() +
plot.males +
plot.stillborn+
plot.spreizer) +
plot_layout(guides = "collect") +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(face = "bold", size = 20),
legend.position = "right",
legend.direction = "vertical")combined
png("./plots/figure-boar.png",
width = 400, height = 300, units = "mm",
pointsize = 10, res = 600)
combined
dev.off()png
2 How to cite R
“All analyses were performed using R Statistical Software (version 4.4.2; R Core Team 2024)”.
Reference: R Core Team (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
citation()To cite R in publications use:
R Core Team (2024). _R: A Language and Environment for Statistical
Computing_. R Foundation for Statistical Computing, Vienna, Austria.
<https://www.R-project.org/>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Manual{,
title = {R: A Language and Environment for Statistical Computing},
author = {{R Core Team}},
organization = {R Foundation for Statistical Computing},
address = {Vienna, Austria},
year = {2024},
url = {https://www.R-project.org/},
}
We have invested a lot of time and effort in creating R, please cite it
when using it for data analysis. See also 'citation("pkgname")' for
citing R packages.version$version.string[1] "R version 4.4.2 (2024-10-31 ucrt)"citation("tidyverse")Um Paket 'tidyverse' in Publikationen zu zitieren, nutzen Sie bitte:
Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R,
Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller
E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V,
Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to
the tidyverse." _Journal of Open Source Software_, *4*(43), 1686.
doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Article{,
title = {Welcome to the {tidyverse}},
author = {Hadley Wickham and Mara Averick and Jennifer Bryan and Winston Chang and Lucy D'Agostino McGowan and Romain François and Garrett Grolemund and Alex Hayes and Lionel Henry and Jim Hester and Max Kuhn and Thomas Lin Pedersen and Evan Miller and Stephan Milton Bache and Kirill Müller and Jeroen Ooms and David Robinson and Dana Paige Seidel and Vitalie Spinu and Kohske Takahashi and Davis Vaughan and Claus Wilke and Kara Woo and Hiroaki Yutani},
year = {2019},
journal = {Journal of Open Source Software},
volume = {4},
number = {43},
pages = {1686},
doi = {10.21105/joss.01686},
}citation("patchwork")Um Paket 'patchwork' in Publikationen zu zitieren, nutzen Sie bitte:
Pedersen T (2024). _patchwork: The Composer of Plots_. R package
version 1.3.0, <https://CRAN.R-project.org/package=patchwork>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Manual{,
title = {patchwork: The Composer of Plots},
author = {Thomas Lin Pedersen},
year = {2024},
note = {R package version 1.3.0},
url = {https://CRAN.R-project.org/package=patchwork},
}citation("lmerTest")To cite lmerTest in publications use:
Kuznetsova A, Brockhoff PB, Christensen RHB (2017). "lmerTest
Package: Tests in Linear Mixed Effects Models." _Journal of
Statistical Software_, *82*(13), 1-26. doi:10.18637/jss.v082.i13
<https://doi.org/10.18637/jss.v082.i13>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Article{,
title = {{lmerTest} Package: Tests in Linear Mixed Effects Models},
author = {Alexandra Kuznetsova and Per B. Brockhoff and Rune H. B. Christensen},
journal = {Journal of Statistical Software},
year = {2017},
volume = {82},
number = {13},
pages = {1--26},
doi = {10.18637/jss.v082.i13},
}citation("car")To cite the car package in publications use:
Fox J, Weisberg S (2019). _An R Companion to Applied Regression_,
Third edition. Sage, Thousand Oaks CA.
<https://www.john-fox.ca/Companion/>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Book{,
title = {An {R} Companion to Applied Regression},
edition = {Third},
author = {John Fox and Sanford Weisberg},
year = {2019},
publisher = {Sage},
address = {Thousand Oaks {CA}},
url = {https://www.john-fox.ca/Companion/},
}citation("emmeans")Um Paket 'emmeans' in Publikationen zu zitieren, nutzen Sie bitte:
Lenth R (2024). _emmeans: Estimated Marginal Means, aka Least-Squares
Means_. R package version 1.10.6,
<https://CRAN.R-project.org/package=emmeans>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Manual{,
title = {emmeans: Estimated Marginal Means, aka Least-Squares Means},
author = {Russell V. Lenth},
year = {2024},
note = {R package version 1.10.6},
url = {https://CRAN.R-project.org/package=emmeans},
}citation("performance")Um Paket 'performance' in Publikationen zu zitieren, nutzen Sie bitte:
Lüdecke et al., (2021). performance: An R Package for Assessment,
Comparison and Testing of Statistical Models. Journal of Open Source
Software, 6(60), 3139. https://doi.org/10.21105/joss.03139
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Article{,
title = {{performance}: An {R} Package for Assessment, Comparison and Testing of Statistical Models},
author = {Daniel Lüdecke and Mattan S. Ben-Shachar and Indrajeet Patil and Philip Waggoner and Dominique Makowski},
year = {2021},
journal = {Journal of Open Source Software},
volume = {6},
number = {60},
pages = {3139},
doi = {10.21105/joss.03139},
}Session Info
sessionInfo()R version 4.4.2 (2024-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)
Matrix products: default
locale:
[1] LC_COLLATE=German_Germany.utf8 LC_CTYPE=German_Germany.utf8
[3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
[5] LC_TIME=German_Germany.utf8
time zone: Europe/Berlin
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] performance_0.13.0 emmeans_1.10.6 car_3.1-3 carData_3.0-5
[5] lmerTest_3.1-3 lme4_1.1-36 Matrix_1.7-1 patchwork_1.3.0
[9] lubridate_1.9.4 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[13] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
[17] ggplot2_3.5.1 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] gtable_0.3.6 bayestestR_0.15.1 xfun_0.50
[4] htmlwidgets_1.6.4 ggrepel_0.9.6 insight_1.0.1
[7] lattice_0.22-6 tzdb_0.4.0 numDeriv_2016.8-1.1
[10] vctrs_0.6.5 tools_4.4.2 Rdpack_2.6.2
[13] generics_0.1.3 datawizard_1.0.0 parallel_4.4.2
[16] pbkrtest_0.5.3 sandwich_3.1-1 pkgconfig_2.0.3
[19] lifecycle_1.0.4 compiler_4.4.2 farver_2.1.2
[22] munsell_0.5.1 codetools_0.2-20 htmltools_0.5.8.1
[25] yaml_2.3.10 Formula_1.2-5 pillar_1.10.1
[28] nloptr_2.1.1 MASS_7.3-61 reformulas_0.4.0
[31] boot_1.3-31 abind_1.4-8 multcomp_1.4-26
[34] nlme_3.1-166 tidyselect_1.2.1 digest_0.6.37
[37] mvtnorm_1.3-3 stringi_1.8.4 labeling_0.4.3
[40] splines_4.4.2 fastmap_1.2.0 grid_4.4.2
[43] colorspace_2.1-1 cli_3.6.3 magrittr_2.0.3
[46] survival_3.7-0 broom_1.0.7 TH.data_1.1-3
[49] withr_3.0.2 backports_1.5.0 scales_1.3.0
[52] timechange_0.3.0 estimability_1.5.1 rmarkdown_2.29
[55] zoo_1.8-12 hms_1.1.3 coda_0.19-4.1
[58] evaluate_1.0.3 knitr_1.49 rbibutils_2.3
[61] mgcv_1.9-1 rlang_1.1.4 Rcpp_1.0.13-1
[64] xtable_1.8-4 glue_1.8.0 see_0.9.0
[67] rstudioapi_0.17.1 minqa_1.2.8 jsonlite_1.8.9
[70] R6_2.6.1